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Abstract 

We explore the impact of obliquity variations on planetary habitability in hypothetical systems with high mutual 
inclination. We show that large-amplitude, high-frequency obliquity oscillations on Earth-like exoplanets can 
suppress the ice-albedo feedback, increasing the outer edge of the habitable zone. We restricted our exploration 
to hypothetical systems consisting of a solar-mass star, an Earth-mass planet at 1 AU, and 1 or 2 larger planets. 
We verified that these systems are stable for 10 8 years with TV-body simulations and calculated the obliquity 
variations induced by the orbital evolution of the Earth-mass planet and a torque from the host star. We ran a 
simplified energy balance model on the terrestrial planet to assess surface temperature and ice coverage on the 
planet’s surface, and we calculated differences in the outer edge of the habitable zone for planets with rapid 
obliquity variations. For each hypothetical system, we calculated the outer edge of habitability for two con- 
ditions: (1) the full evolution of the planetary spin and orbit and (2) the eccentricity and obliquity fixed at their 
average values. We recovered previous results that higher values of fixed obliquity and eccentricity expand the 
habitable zone, but we also found that obliquity oscillations further expand habitable orbits in all cases. 
Terrestrial planets near the outer edge of the habitable zone may be more likely to support life in systems that 
induce rapid obliquity oscillations as opposed to fixed-spin planets. Such planets may be the easiest to directly 
characterize with space-borne telescopes. Key Words: Exoplanets — Habitable zone — Energy balance models. 
Astrobiology 14, 277-291. 


1. Introduction 

T he habitability of a world depends on a host of prop- 
erties, from observable quantities like its mass and dis- 
tance from the parent star to those that are difficult if not 
impossible to measure: atmospheric composition, surface 
reflectivity, ice, water distribution, and so on. In the case of 
stars as massive as our Sun, detecting Earth-mass planets in 
any orbit is difficult with modern technology. In the last de- 
cade, attention has turned primarily to the discovery of rocky 
planets orbiting in the habitable zone (HZ), a shell around a 
luminous object in which an Earth-like planet could support 
liquid water on its surface (Dole, 1964; Kasting et al., 1993; 
Kopparapu et al., 2013), as these worlds are best suited for the 
development and sustainment of life as we know it. In its 
latest revision (Kopparapu et al., 2013), the HZ is calculated 
for a highly idealized case, in which many properties of the 
star, planet, and planetary system are ignored. Following the 
identification of possible processes that can impact habit- 
ability (Heller and Armstrong, 2014), we explore how grav- 


itational perturbations from additional planets can affect the 
climate. We find that, in many cases, these perturbations can 
extend the outer edge of the HZ, thereby increasing the 
number of planets in the Galaxy that are potentially habitable. 

While the vast majority of work on habitability has used a 
replica of Earth to determine orbits that are potentially 
habitable, there are notable exceptions. Some studies have 
explored the habitability of synchronously rotating planets 
(Joshi et al., 1997; Joshi, 2003; Edson et al., 2011; Pierre- 
humbert, 2011; Wordsworth et al., 2011; Yang et al., 2013). 
Others (Abe et al., 2011; Zsom et al., 2013) considered 
planets that are much drier than Earth. Several studies 
(Williams and Pollard, 2002, 2003; Dressing et al., 2010; 
Spiegel et al., 2010) varied the eccentricity and obliquity of 
an Earth-like planet and found that larger values tend to 
increase the globally averaged temperature on a planet, 
while holding the semimajor axis constant. 

While these studies made great strides in understanding 
Earth’s climate sensitivity to rotation rate, obliquity, and 
eccentricity, aside from Spiegel et al. (2010), they largely 
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ignored that the latter two properties evolve with time due to 
gravitational perturbations from other bodies. Earth main- 
tains a relatively constant axial tilt either due to the presence 
of the Moon, as suggested by Laskar et al. (1993), or due to 
the inherent stability of Earth’s axis, as indicated by Lis- 
sauer et al. (2012). However, as illustrated by Williams 
(1998a, 1998b), changes in the architecture of our solar 
system — such as moving Jupiter inward — can result in 
dramatic variations in the obliquity of Earth even with the 
presence of a Moon. Still, it has been suggested that the 
relatively small variations in Earth’s obliquity result in a 
stable climate conducive to the development of life. Adding 
to this stability is the fact that the orbital eccentricity re- 
mains smaller than about 0.05 due to the approximately 
circular orbits of the large planets of the Solar System. 

It is possible that small obliquities and circular orbits are 
not a requirement for habitability. Williams and Pollard 
(2003) used general circulation models (GCMs) to determine 
how different obliquity variations affect Earth’s climate. 
They found that Earth-like planets with high obliquities were 
no more likely to experience extreme runaway greenhouse 
or snowball Earth events, making them just as habitable as 
Earth. Later, Spiegel et al (2010) examined how large ec- 
centricity oscillations affect the climates of rocky exoplanets. 
They found that, in some cases, planets could break out of a 
snowball event during periods of high eccentricity. It was the 
goal of the current study to build on these previous results and 
explore self-consistent models of the climates of planets that 
experience rapid, large-amplitude, and possibly chaotic os- 
cillations of eccentricity and obliquity. 

Orbit-induced seasonal effects like the ice-albedo feedback 
determine the limit of the outer edge of the HZ. As the surface 
temperature drops, volatile ices such as C0 2 and water can 
condense on the surface. The high albedos of their solid 
phases inhibit a planet’s ability to absorb solar radiation, 
which reduces the temperature further. Mars, if it possessed 
sufficient surface water, would have been in danger of falling 
prey to these snowball episodes. Geological evidence exists 
for these episodes in Earth’s past (Hoffman et al., 1998). For 
Earth, a dynamic C0 2 recycling system works to offset the 
negative effects of these events on timescales of millions of 
years (Walker et al., 1981). As the planet cools, weathering 
rates slow down and lock C0 2 in the atmosphere. As the C0 2 
builds up, the greenhouse effect increases, eventually melting 
the ice. On early Mars, as now, such robust C0 2 cycling could 
not have been enough to resurrect the planet from these 
snowball events. From seasonally resolved modeling, it seems 
that other factors, in particular orbital and obliquity variations, 
would have to play a role (Armstrong et al., 2004). Mars’ 
obliquity has probably undergone significant evolution in the 
past (Laskar and Robutel, 1993) due to gravitational torques 
by the other planets in the Solar System, a property that may 
have permitted at least episodic liquid water at the surface. 

With this Solar System context in mind, we turn our at- 
tention to predicting the habitability of exoplanets. We do 
not anticipate the photometric and spectroscopic data nee- 
ded to characterize exoplanets until the launch of the James 
Webb Space Telescope in ~2017. Spacecraft like Kepler 
and ground-based radial velocity surveys are already re- 
turning data that can pin down, through computational 
analysis, important orbital parameters that impact the cli- 
mate, such as eccentricity, timing of perihelion passage, and 


the evolution of the spin axis of the planet. Through the 
coupling of these data to A-body simulations and simple, 
fast climate codes, a more comprehensive picture of the HZ 
of a system can be obtained. This type of analysis can be 
used to prioritize future characterization observations, which 
are likely to be challenging, expensive, and based on pre- 
cious little information. In the present study, we found that 
planetary system architecture, that is, the distribution of 
masses and orbits of the other planets in the system, can play 
a significant role in defining the extent of the HZ. 

The rotational evolution of the bodies in our solar system 
can be accurately modeled because the masses and orbits of 
the planets are extremely well measured. For exoplanets, the 
situation is more difficult. Radial velocity surveys ( e.g ., 
Butler et al., 2006) are only able to place a lower bound on 
mass and cannot measure the relative inclination between 
orbital planes. Kepler can constrain inclinations (Fabrycky 
et al., 2012) but is heavily biased toward the discovery of 
planets in coplanar configurations. While these systems are 
analogous to our solar system, we note that for the one 
system for which astrometry (which is not biased toward 
any particular inclination) has measured a mutual inclina- 
tion, v Andromedae (McArthur et al., 2010), the relative 
inclination is 30 degrees. Moreover, studies that predict the 
large eccentricities of exoplanetary orbits simultaneously 
predict large inclinations (Marzari and Weidenschilling, 
2002; Chatterjee et al., 2008; Raymond et al., 2010; Barnes 
et al., 2011). It is possible, perhaps likely, that there exists a 
population of planetary systems with large inclinations and a 
potentially habitable planet. These architectures will induce 
much larger changes in orbital inclination, which in turn 
induces large obliquity oscillations. The Gaia mission may 
be able to determine the range of architectures for giant 
planets (Casertano et al., 2008; Sozzetti et al., 2013). 

As no rocky planet is currently known to orbit with sib- 
ling planets with high mutual inclinations, we explored the 
phenomena with 17 hypothetical, dynamically stable sys- 
tems. We found that there is a direct link between the orbital 
architecture of a planetary system and the possible range 
of climate conditions on a potentially habitable planet. We 
used these models to constrain the orbital conditions of a 
hypothetical planet and found that orbital and rotational 
evolution tends to push the outer edge of the HZ out, relative 
to planets where no evolution occurs. 

Below, we outline a model that links physically realistic 
orbital architectures to the spin evolution of a hypothetical 
Earth-like planet and finally to its climate. In Section 2, we 
discuss the motivation behind the systems we have modeled 
in an effort to obtain a set that spans a range of orbital 
elements. In Section 3, we outline the model used to evolve 
the spin axis of the planet. In Section 4, we present a 
simplified energy balance model (EBM) designed to be 
robust across wide variations of these orbital changes and 
fast enough that million-year integrations require only a few 
hundred seconds of computational time. We then present 
the results of these models in Section 5 followed by a dis- 
cussion in Section 6. 

2. Orbital Simulations 

To make a first assessment of the potential of similar 
architectures to support habitable worlds, we created 17 
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Table 1. Summary of 17 Systems Modeled 
Including the Mean Semimajor Axis, a ; 
Mean Eccentricity, e; Mean Inclination, i; 
and Planetary Mass, M p , in Earth Masses 


System 

N 

a, AU 

e 

i, deg 

M, M® 

1 

1 

1.0 

0.05 

1.4 

0.3 


2 

5.2 

0.05 

0.9 

332.9 


3 

9.5 

0.05 

1.7 

133.2 

2 

1 

1.0 

0.14 

17.5 

1.0 


2 

0.4 

0.30 

19.2 

10.0 


3 

3.0 

0.30 

7.9 

10.0 

3 

1 

1.0 

0.33 

20.2 

1.0 


2 

10.0 

0.08 

17.7 

3178.0 


3 

19.5 

0.07 

12.5 

3178.0 

4 

1 

1.0 

0.09 

11.3 

1.0 


2 

0.1 

0.30 

10.9 

3178.0 


3 

34.0 

0.64 

2.2 

3178.0 

5 

1 

1.0 

0.08 

30.1 

1.0 


2 

5.0 

0.30 

0.0 

3178.0 

6 

1 

1.0 

0.10 

10.0 

1.0 


2 

5.0 

0.30 

0.0 

317.8 

7 

1 

1.0 

0.08 

30.0 

1.0 


2 

5.0 

0.30 

0.0 

317.8 

8 

1 

1.0 

0.0001 

15.3 

1.0 


2 

0.4 

0.00002 

19.6 

10.0 


3 

2.5 

0.0001 

8.6 

10.0 

9 

1 

1.0 

0.02 

2.3 

1.0 


2 

14.9 

0.07 

11.7 

3178.0 


3 

29.3 

0.09 

8.3 

3178.0 

10 

1 

1.0 

0.08 

21.2 

1.0 


2 

0.1 

0.30 

20.5 

3178.0 


3 

34.2 

0.64 

3.9 

3178.0 

11 

1 

1.0 

0.08 

32.7 

1.0 


2 

0.1 

0.30 

30.8 

3178.0 


3 

34.7 

0.64 

5.6 

3178.0 

12 

1 

1.0 

0.19 

20.8 

1.0 


2 

0.6 

0.21 

21.4 

1.0 


3 

2.5 

0.30 

23.3 

1.0 

13 

1 

1.0 

0.02 

7.2 

1.0 


2 

10.0 

0.06 

6.4 

317.8 


3 

29.9 

0.04 

3.7 

317.8 

14 

1 

1.0 

0.02 

14.3 

1.0 


2 

10.0 

0.06 

12.7 

317.8 


3 

29.9 

0.04 

7.3 

317.8 

15 

1 

1.0 

0.19 

21.9 

1.0 


2 

10.0 

0.06 

19.1 

317.8 


3 

29.9 

0.04 

10.9 

317.8 

16 

1 

1.0 

0.02 

14.0 

1.0 


2 

5.0 

0.04 

10.8 

317.8 


3 

10.0 

0.06 

19.3 

127.1 

17 

1 

1.0 

0.08 

30.0 

1.0 


2 

5.0 

0.30 

0.1 

127.1 


The systems consist of either two or three planets, one being an 
Earth-mass planet located at 1.0 AU. 


hypothetical systems with moderate inclinations, always 
initially including an Earth-like planet on a circular orbit 1 
AU from a solar-mass star. The orbital architectures are 
arbitrary but consistent with the distribution of orbital ele- 
ments of known planets. While the potentially habitable 
world is always the same, its siblings have a wide array of 
properties. Systems consist of 2-3 planets, with eccentrici- 
ties ranging from 0 to 0.3 and mutual inclinations from 10 to 


30 degrees. The orbital properties of these cases, in astro- 
centric coordinates, are presented in Table 1. 

These systems were selected after careful consideration 
of orbital stability. We numerically integrated each case 
for 100 Myr in “hybrid” mode with Mercury (Chambers, 
1999) and confirmed that the evolution of every orbital 
element appeared periodic. Additionally, energy was con- 
served to better than 1 part in 10 6 , which is sufficient for 
numerical accuracy (Barnes and Quinn, 2004). We then 
rotated each system such that the reference plane corre- 
sponds to the fundamental plane. For those systems with 
super-Jupiter planets, this conversion can change some 
orbital elements significantly. We then reran each case at 
very high resolution for ~ 1 Myr, conserving orbital an- 
gular momentum to 1 part in 10 12 . We are therefore con- 
fident that no numerical inaccuracies are propagated into 
the rotational calculations. 

3. Obliquity Modeling 

Using the results from the orbital runs described above, 
we employed the obliquity model of Laskar (1986a,b) as 
used in previous orbit coupled modeling by Armstrong 
et al. (2004). There are two primary factors that influence 
the evolution of the obliquity. First, variations in the ge- 
ometry of the orbit that are governed by the overall system 
architecture are present in the orbital elements derived 
from the A-body simulations. The gravitational influence 
of the other massive bodies in the system affects the eccen- 
tricity, e, the inclination, i, the argument of perihelion, co, 
and the longitude of the acceding node, £?, of the Earth-sized 
planet, here cast in terms of the eccentricity-inclination 
variables 


h = e sin (m) 

a) 

k = e cos (m) 

(2) 

' ,=sin (£ 

) sin (( 2 ) 

(3) 

«= sin ({) 

1 cos (( 2 ) 

(4) 


where w = Q + co is longitude of perihelion. Changes in 
parameters like the inclination result in changes in the spin 
axis orientation relative to the fundamental plane, that is, the 
obliquity, as illustrated in Fig. 1. In other words, the rota- 
tional angular momentum is decoupled from the orbital 
angular momentum. 

In addition to these geometric factors, the direct torques 
from the central star are included as a term, R(ij /), in the 
precession, p A , and obliquity, i/a evolution equations, 

-^7 = R( i/0 - cot (i/0 [A(p, q ) sin (p A ) 
at (5) 

+ B(p, q) cos ( p A )] -2 C(p,q) -p g 


— = - B(p, q) sin ( p A ) + A(p, q) cos ( p A ) ( 6 ) 

at 
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FIG. 1 . A simplified schematic illustrating 
how the evolution of the inclination leads to 
an evolution in obliquity. As the inclination 
of the orbit increases, the spin axis continues 
to point at a fixed position in space, causing 
the angle between the spin axis and the fun- 
damental plane to increase. In this figure, the 
inclination changes in such a way that the 
obliquity goes from a starting value of 23.5 
degrees to 0 degrees. (Color graphics avail- 
able online at www.liebertonline.com/ast) 


with 


3 k 2 M 

R( ip)=— — E d S 0 cos(i/0 

7 / 


(7) 


A(p, q) = 


\J\ -p 1 -q 2 


[q-pC(p,q )] (8) 


B(p, q) = 


\/ l 


= [p-qC(p,q )] (9) 


C(p,q) = qp-pq (10) 


So= 1(1 —e 2 ) - 3/2 — 0.522 x 10“ 6 (11) 

Here, the semimajor axis, a , is measured in AU; the planet’s 
angular velocity, v, is measured in rad day _ M= 1.0 in solar 
units; k 2 = GM/4n 2 , with G, the gravitational constant, mea- 
sured in units of AU 3 M~ l day 2 ; and E D is the dynamical 
ellipticity, a measure of the non- sphericity of the planet. Fi- 
nally, the relativistic precession is accounted for with 


k r 

Pg= 2(1 -e 2 ) 
where the value of K r is 


( 12 ) 


namical ellipticity and orbits a Sun-like star. The obliquity, 
precession angle, and rates of change in these parameters are 
computed each time step from the input orbital parameters, 
and those results are used as the initial conditions for the 
following time step in which a fourth-order Runge-Kutta in- 
tegrator is used. These results are then available for further 
coupling to the EBM described below. The primary limitation 
to the method is that only direct torques from the central body 
are included. Also, there is an implicit assumption that the 
planet is rapidly rotating and that it can be accurately described 
by the dynamical ellipticity relative to the central body. Since 
none of our models involve tidally locked worlds or close- 
passing planets or moons, these assumptions are justified. 

Of the 17 systems studied, we selected seven to analyze 
that represented the full spectrum of outcomes from the 
models: 

(1) Earth- Jupiter- Saturn system for comparison (System 
1, Baseline — Fig. 2). 

(2) Two systems with high mean obliquity but modest 
variations (Systems 2 and 3 — Figs. 3 and 4). 

(3) Two systems with wide and rapid variations in 
obliquity (Systems 4 and 5 — Figs. 5 and 6). 

(4) Two systems with wide and slow variations in 
obliquity (Systems 6 and 7 — Figs. 7 and 8). 

Figure 2 illustrates the baseline run with (Moon-less) 
Earth, Jupiter, and Saturn. The panels on the left represent 
the relevant orbital parameters used in the obliquity calcu- 
lations, and the right-hand panels are the parameters of in- 
terest for the climate calculations. These solutions are 
qualitatively similar to recent work ( e.g ., Lis sauer et al., 


n a 

kr= c 2 (l+M p /M) (13) 

where c is the speed of light, M p is the planet’s mass in units 
of solar masses, and n and a are related via Kepler’s law 

n 2 a 3 = k 2 (l +M p /M) (14) 

as outlined in Laskar (1986a, 1986b). The values used for 
these model parameters are listed in Table 2. 

For each of the 17 orbital runs, the Earth-mass planet in 
the simulation has a 1-day rotation with an Earth-like dy- 


Table 2. Parameters for the Obliquity Calculation 


Parameter 

Symbol 

Value 

Units 

Gravity constant 

k 

0.01720209895 

AU 3 Mq 1 day “ 2 

Angular spin 

V 

2n 

rad day -1 

rate 

Dynamical 

E D 

0.00328005 

unitless 

ellipticity 

Initial obliquity 


23.44 

degrees 

Initial 

Pa 

0.0 

degrees 

precession 

Planet mass 

M p 

1.0 

Earth masses 
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FIG. 2. Orbital-rotational results 
for System 1, Baseline, the Earth- 
like comparison system. The left 
column shows the variations of ec- 
centricity, inclination, and longi- 
tude of ascending node; the right 
column shows the variations of the 
argument of perihelion, obliquity, 
and precession rates. 



2012) but differ in the detailed periodicities and magnitudes 
of the obliquity and precession due to the specific orbital 
solution used in the model. However, the climatically im- 
portant parameters have relatively low amplitude and are 
slowly varying compared to other models. 

Figure 3 shows a model with the Earth-mass planet in a 
system with two other planets 10 times its size, with large 
eccentricity variations and high mutual inclination. This 
architecture results in wide swings in obliquity from the 
starting value to about 80 degrees but has relatively slow 
variations. 


Figure 4 shows a model with two 10 Jupiter-mass planets 
orbiting in the system with the Earth-mass planet, again with 
high eccentricity, high mutual inclinations, and an obliquity 
that oscillates around 85 degrees with a period of approxi- 
mately 200,000 years. 

Figures 5 and 6 illustrate some of the cases where the 
obliquity variations become rapid (periods of < 16,000 
years) at relatively high amplitudes between 10 and 60 de- 
grees. In each of these cases, there is high mutual inclination 
and large eccentricities due to two 10 Jupiter-mass planets in 
Fig. 4 and one Jupiter-mass planet in Fig. 5. 


FIG. 3. Orbital-rotational results 
for System 2. The left column 
shows the variations of eccentricity, 
inclination, and longitude of as- 
cending node; the right column 
shows the variations of the argu- 
ment of perihelion, obliquity, and 
precession rates. 
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FIG. 4. Orbital-rotational results 
for System 3. The left column 
shows the variations of eccentricity, 
inclination, and longitude of as- 
cending node; the right column 
shows the variations of the argu- 
ment of perihelion, obliquity, and 
precession rates. 




Lastly, Figs. 7 and 8 show an Earth-like planet that has a 
Jupiter-mass planet in the system with a mutual inclination 
of 10 degrees (Fig. 7) and 30 degrees (Fig. 8), resulting in a 
wide range in obliquity that varies more slowly than the 
previous cases. 

The precession rates illustrated in the figures are positive, 
as expected, except near {// = 0 or ^>90 degrees. These are 
numerical artifacts of the geometry of the semianalytical 
model. As the obliquity goes beyond 90 degrees, the pre- 
cession rate reverses due to the fact that the “south” pole is 
now the “north” pole. Additionally, the cotangent term in 


Eq. 5 makes the precession poorly defined near 0 degrees. 
From a climate perspective, the precession is unimportant at 
low obliquities, and this is only relevant for “bookkeeping.” 

4. A Simplified Energy Balance Model 

With the orbital variations and obliquity calculations in 
hand, we assessed the surface conditions using a simplified 
EBM. This calculation requires a model that can simulate 
planets on the timescales for glacier growth/retreat: 10 4 to 
10 5 years. GCMs are generally too complex to be practical 



400 



0.0 0,2 0.4 0.6 0.8 1.0 



IS I 1 1 1 1 

0.0 0,2 0,4 0,6 0,8 1.0 



0,0 0-2 0.4 0.6 0-8 10 


FIG. 5. Orbital-rotational results 
for System 4. The left column 
shows the variations of eccentricity, 
inclination, and longitude of as- 
cending node; the right column 
shows the variations of the argu- 
ment of perihelion, obliquity, and 
precession rates. 
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FIG. 6. Orbital-rotational results 
for System 5. The left column 
shows the variations of eccentricity, 
inclination, and longitude of as- 
cending node; the right column 
shows the variations of the argu- 
ment of perihelion, obliquity, and 
precession rates. 


30.30 


30,25 

30,20 

30.15 

30.10 


- 30.05 
30.00 
29,95 

29.90 

0 






on these timescales. And both GCMs and detailed EBMs, to 
a greater or lesser extent, require detailed information about 
the planet or assumptions based on relevant observations. In 
the case of exoplanets, especially those with widely varying 
orbital parameters, making these assumptions is difficult if 
not impossible. To include enough detailed physics to allow 
for direct calculation of climate effects in absence of ob- 
servation would render our code too cumbersome. To act as 
a first-order comparison between the other systems and the 


baseline model, with as few input parameters as possible, we 
have developed a simplified EBM to examine the conditions 
on the surface. While the lack of detail makes specific cli- 
mate predictions difficult for individual planets, general 
comparisons should be valid. 

We modeled the surface as a one-dimensional latitude 
grid of 90 bands from - 90 S to + 90 N and the atmosphere 
as a single slab gray absorbing layer. We modified a sim- 
plified model outlined by Armstrong et al. (2004) to include 




FIG. 7. Orbital-rotational results 
for System 6. The left column 
shows the variations of eccentricity, 
inclination, and longitude of as- 
cending node; the right column 
shows the variations of the argu- 
ment of perihelion, obliquity, and 
precession rates. 
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FIG. 8. Orbital-rotational results 
for System 7. The left column 
shows the variations of eccentricity, 
inclination, and longitude of as- 
cending node; the right column 
shows the variations of the argu- 
ment of perihelion, obliquity, and 
precession rates. 


basic deposition and evaporation of water ice over the 
seasonal cycle. 

For each time step in the million-year orbital-obliquity 
model, we modeled three years of the seasonally resolved 
climate, one year to act as a spin-up from the initial con- 
ditions, and two years to compute global surface properties. 
Each model year starts with t = 0 at the time of perihelion 
passage, r peri = 0. From this, we can compute the eccentric 
anomaly, E , from the implicit equation 

( c? y /2 

t ^=(gmJ (E ~ esinE) (15) 

using Newton’s method (where all the values above are 
measured in MKS units). The eccentric anomaly is related to 
the true anomaly, /, by 

l+eV /2 


Id = — [77 sin <5* sin S + sin (77) cos S * cos (5] (19) 
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where 77 is the half- angle of daylight — a measure of the 
length of the day — at a given latitude 5; and the substellar 
latitude, <5* = ^ sin(L s ), is determined by the obliquity, \//, 
and the stellar longitude, L s = L sp +/ where L s = 0 is the 
northern hemisphere vernal equinox. L sp = L sp0 - p A - co is 
the solar longitude of perihelion, determined by the arbitrary 
constant L sp0 (in our case, 90 degrees represents northern 
hemisphere summer solstice), the spin precession, p A , and 
argument of perihelion, co. The half-angle of daylight, r\ , is 
given by 

cos 77 = — tan S tan (5* \5\ < 90— |<5*| (20) 


cos 77= -1 S< —90 + <5* or £<90 + <5* (21) 


/ = 2 arctan 



With the true anomaly, we compute the instantaneous dis- 
tance to the star, 


a( 1 — e 2 ) 

1 + e cos f 


(17) 


The instantaneous stellar distance allows us to compute the 
incident stellar radiation, 


S p = 


4nr 2 


(18) 


where L* is the stellar luminosity. From this, we can com- 
pute the daily mean top-of- atmosphere instellation at any 
point on the globe, 


cos 77=1 <5 >90 — <5* or £<90 - <5* (22) 

To estimate the surface temperature, we first determine the 
energy balance at the surface at each latitude bin, 

AE = I d (\ -A)- e s (T r e 4 + F surf (23) 

where A is the surface albedo of either ground or ice 
depending on local conditions, s s is the surface emis- 
sivity, T e is the planet’s atmosphere-free equilibrium 
temperature, and F surf is the heat flux from the surface. 
The atmosphere is modeled as a single slab with opacity 
t that is equivalent to the number of absorbing layers 
required to achieve a surface temperature of T s , related 
to T e by 
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r s = (i + T) 1/4 r e (24) 

The t parameter is essentially a one-dimensional model of the 
greenhouse effect. That is, adjusting t as a free parameter 
allows us to include a contribution to greenhouse warming, 
with t = 0.095 reproducing the Earth-like global surface 
temperature in our baseline model. To estimate the surface 
temperature, we set AE=0, solve first for the equilibrium 
temperature, and then solve for the surface temperature. 

Ice deposition is handled parametrically in the model by 
choosing a global ice deposition rate that comes into effect 
when the surface temperatures dip below 273 K. At this 
point, the surface albedo is set to the value for ice, and ice is 
allowed to accumulate. If the surface temperature exceeds 
273 K in the presence of ice, the surface temperature is held 
at 273 while the deposits are evaporated according to the 
difference between what the temperature of the surface 
would be in the absence of ice and the freezing point of 
water, T ice = 273. In this case, the mass loss rate in units of 
kg s' 1 m -2 is given by 

At L h 1 ’ 

where L h = 3.34x 10 5 is the latent heat of fusion for ice in 
J kg -1 . Once the ice is gone, the albedo is reset to the 
surface value, and the surface temperature evolves normally. 
The ice deposition rate is set such that the thickness and 
extent of the ice- snow regions are roughly comparable with 
values for present-day Earth in the baseline model. 

In summary, the only “parameters” are the albedo of 
land/ice, the atmospheric opacity (which is varied as a free 
parameter in the study), the ice deposition rate, and a pa- 
rameterized version of thermal exchange with the surface. 
To calibrate the model, we chose reasonable values for the 
albedo, tuned the deposition rate to get reasonable ice caps, 
and modified the atmospheric opacity and surface flux to hit 
a global mean temperature of ~288 K. The parameters of 
the model are outlined in Table 3. 

This approach has several important features: 

(1) The model is extremely fast, allowing us to run 
thousands of cases over millions of years to explore 
the full parameter space of the model. 


Table 3. Parameters for the Climate Calculation 


Parameter 

Symbol 

Value 

Units 

Stellar mass 

M* 

1.0 

Solar masses 

Stellar luminosity 

L* 

1.0 

Solar luminosities 

Land albedo 

^land 

0.4 

unitless 

Ice albedo 

A. 

0.6 

unitless 

Surface heat flux 

^surf 

88 

watts 

Baseline opacity 

T 

0.095 

unitless 

Surface emissivity 

^surf 

1.0 

unitless 
kg m -2 s' 1 

Ice/snow 
deposition rate 

^ice 

5. Ox 10 -5 


The parameters of the climate model are selected to be within 
reasonable limits and still reproduce climates similar to that of Earth 
under current conditions. The only other free parameters are the 
atmospheric opacity and the distance from the star, both varied as 
part of the study. 


(2) The physics is extremely simplified and coupled 
tightly to the orbital parameters. Therefore, we can 
see the first-order effects of the orbital variations, the 
main goal of this study. 

(3) The model is determined by only three free parame- 
ters: the opacity (or number of absorbing layers), the 
surface heat flux, and the ice deposition rate. 

(4) The final results are normalized to the baseline run, 
which, while limiting what we can say about the 
specific climate, allows for a quantitative comparison 
of the outer edge of the HZ for the hypothetical 
planets and modern Earth. 

This model reproduces the first-order climate effects and 
the ice-albedo feedback (see Fig. 9). Each panel plots the 
two model years after the spin-up year. The top plot is the 
baseline run, showing the surface temperatures for our 
nominal “Earth.” The discontinuities result from the abrupt 
change in surface albedo. The second model is identical to 
the first, except the eccentricity has been increased to 0.15 to 
show the asymmetric effect on the climate, in this case 
creating a more persistent polar cap in the southern hemi- 
sphere. The third model is the same as the first but shows the 
effect of increasing the obliquity to 90 degrees. 

Using this model, we computed approximately 4400 test 
cases for the 17 models over a range of semimajor axes from 
0.80 to 3.0 and opacities ranging from 0.095, obtained from 
tuning the Earth baseline model to 288 K, to 0.26, nearly 3 
times the amount of absorbing material in the atmosphere. 
For each time step in the orbital model, we ran an instance 
of the EBM, computing the three-year climate run, and 
averaged the last two years to compute the global mean 
temperature. In addition, we computed two additional 
measures of the habitability of the surface, the Temperature 
Habitability Index (THI), and the Ice Habitability Index 
(IHI). THI is defined here as the time-averaged fraction of 
the surface that is between 273 and 373 K. IHI is defined as 
the time-averaged ice-free regions of the planet regardless of 
surface temperature. We then normalized the quantities to 
the THI and IHI for the t- 0 baseline Earth, which has a 
THI = 0.766 and IHI = 0.485. If the THI (or IHI) is greater 
than 1.0, this planet is “more habitable” than the baseline 
Earth. If the THI or IHI is less than 1 .0, it is less habitable. A 
THI of 0 indicates a world that is either frozen (in which 
case the IHI will also be 0) or has no temperatures below 
373 K. The maximum values for THI and IHI by this metric 
are 1.3 and 2.1, respectively. 

This method has some distinct advantages for exploring 
the effects of the orbit on climate. Since the model is sea- 
sonally resolved, even some planets that might have an 
annually averaged global mean surface temperature below 
freezing may still experience significant periods of time in 
the “habitable” range and, due to the orbital effects, avoid a 
planetwide snowball. 

5. Results 

Figure 10 shows the aggregate results for the seven sys- 
tems studied in detail. The first column is the THI, the 
second column is the IHI, and the third column is the global 
mean temperature, all as a function of both the distance from 
the host star and the opacity of the atmosphere for Systems 
1-7. Since our simplified EBM ignores complications due to 
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FIG. 9. Baseline climate models for an Earth-like planet 
with i/4 = 25, e = 0.0 (top); 1/4 = 25, e = 0.15 (middle); and if/ = 90, 
<? = 0.0 (bottom). The discontinuities are caused by the change 
in albedo between an ice/snow-covered and ice-free surface. 


the runaway greenhouse, we will restrict our discussion to 
the outer edge of the HZ. 

The baseline model shows that the outer edge of the HZ is 
approximately 1.4 AU for Earth at t = 0.095. Each of the 
other high-obliquity systems and/or high-eccentricity systems 
shows a systematic — and sometimes dramatic — increase in 
the outer edge of the HZ. Those increases are tabulated in 
Table 4. The maximum increase is System 3, which has both 
high eccentricity and large obliquity, with large variations in 
both parameters. 

To separate the effects due to the large amplitudes of the 
eccentricity and obliquity and their variations, we performed 
static calculations using the mean values of e and 1/4 to 
compute the static outer edge. By comparing this value to 
the outer edge of the HZ in the variable runs, / out , we can 
determine how much of the expansion of the outer edge is 
due to the variability and how much is due to the large 
values of e and 1/4. 

Table 4 lists this comparison as the HZ enhancement 
factor, as a percentage, due to the static orbital properties of 
the simulations, 


- lj x 100 (26) 

and the HZ enhancement factor due to the variability of 
those parameters, 


Es 




Vbaseline 


Ey = (j — — l) x 100 (27) 

\ /static J 

where /baseline is the outer edge of the baseline system, and 
/static is the outer edge in the static case. From this analysis, 
we see that the increase in the outer edge is dominated by 
the large values of e and 1/4. However, a non-negligible 
component — in one case the sole component, as in System 
7 — is due to the variability of those values. 

The values for E s and E w are listed in Table 4 and shown 
schematically in Fig. 1 1. In the figure, the height of the bar is 
equal to E s for each system. The green portion of the bar 
indicates the percentage of the enhancement that is due to the 
variability alone. For example, our most enhanced system, 
System 3, had a 93% increase in the HZ compared to the 
baseline model. Of that, 8% of the enhancement is due to the 
variability of the system. In one case, System 7, the en- 
hancement is entirely due to the variability of the parameters. 

In some cases, the nonvariable systems move the HZ 
inward from the baseline value of 1.4. System 9 shows no 
enhancement in the full simulation, but the HZ moves in- 
ward when variability is removed (hence the negative 
“enhancement”). In System 17, the static run produces an 
outer limit that is 7% smaller than the baseline system 
(again causing negative enhancement), but dynamics allow 
a 15% increase compared to its static value. 

Figure 12 compares the outer edge of the HZ for the vari- 
able cases (top panels) and static cases (bottom panels) as a 
function of the eccentricity and obliquity. The error bars on 
Systems 1-17 are derived from the standard deviation of the 
eccentricity and obliquity from the complete simulations. The 
points in the top panels lie to right of those in the bottom, 
showing that the variability increases the HZ. However, the 
effect is most strongly correlated with the obliquity. 
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FIG. 10. The temperature habitability index (THI, left column), the ice habitability index (IHI, middle column), and the 
mean global temperature (right column) for seven systems listed in Section 3. From top to bottom: System 1 (Baseline) 
through System 7. 


In an effort to quantify the relationship, Table 5 lists the 
linear regression slope and intercept, along with an error- 
weighted goodness of fit indicator, 

X 2 =Z— ( 28 ) 

where y t is the “observed” value of either i/y or e,/(xO is the 
linear fit, and is the uncertainty computed for the given 
parameter. In the nonvariable cases, since we have no es- 
timate of the uncertainty, the errors in the goodness of fit 


calculation are taken as 1% of the value of the data point. 
The results show little correlation with eccentricity but a 
very strong relationship between obliquity and the outer 
edge of the HZ. Removing the variability reduces the in- 
tercept by 5 degrees and steepens the slope by 7 degrees per 
astronomical unit, which means the outer edge systemati- 
cally moves inward when the variability is removed. 

6. Discussion 

Our simulations show that the evolution of planets’ orbit 
and rotation can increase the maximum separation between 
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Table 4. Summary of the Systems Used in the Complete Climate Comparison 


N 

e 

P e , Myr 

•A, deg 

P,l„ Myr 

1 out > AU 

1 static > AU 

E s , % 

E v , % 

1 

0.05 + 0.023 

0.25 

22.6 ± 1.1 

0.17 

1.4 

1.4 

n/a 

n/a 

2 

0.14 ±0.042 

0.20 

60.6 ±20.2 

0.33 

2.1 

2.0 

50 

5 

3 

0.33 + 0.025 

0.04 

78. 6± 18.7 

0.17 

2.7 

2.5 

93 

8 

4 

0.08 ±0.036 

0.02 

40.7 ± 12.1 

0.02 

1.9 

1.7 

36 

12 

5 

0.08 ±0.038 

0.05 

40. 6± 19.1 

0.02 

1.9 

1.7 

36 

12 

6 

0.10 ±0.048 

0.14 

58.7±21.4 

0.33 

2.1 

1.9 

50 

11 

7 

0.08 ±0.041 

0.07 

28. 8± 11.4 

0.17 

1.7 

1.4 

21 

21 

8 

0.0001 ±0.00006 

0.007 

37.6± 8.6 

0.23 

1.7 

1.6 

21 

6 

9 

0.09 ±0.038 

0.01 

1 8.4 ± 4. 1 

0.02 

1.4 

1.3 

-7 

8 

10 

0.02 ±0.001 

0.08 

23.9 ±0.7 

0.04 

1.4 

1.4 

0 

0 

11 

0.08 ±0.035 

0.01 

34.3 ±7.1 

0.02 

1.7 

1.6 

21 

6 

12 

0.19 ±0.064 

0.33 

34.3 ±7.4 

0.20 

1.9 

1.6 

36 

19 

13 

0.03 ±0.010 

0.33 

23.7±0.6 

0.09 

1.4 

1.4 

0 

0 

14 

0.02 ±0.008 

0.25 

24.2 ±1.0 

0.09 

1.4 

1.4 

0 

0 

15 

0.17 ±0.026 

0.25 

22.6 ±1.2 

0.08 

1.5 

1.4 

7 

7 

16 

0.02 ±0.008 

0.33 

47. 9± 19.8 

0.20 

1.9 

1.8 

36 

6 

17 

0.08 ±0.043 

0.17 

24.6 ±4.7 

0.10 

1.5 

1.3 

-7 

15 


Table includes the system number, N; the mean eccentricity, e, along with its dominate period of oscillation, P e , in millions of years; 
mean obliquity, i/q with its period, P x j J ; the calculated outer edge of the habitable zone, / out , for temperature based on the habitability index; 
the outer edge of the habitable zone for the static runs, / static ; and the orbit enhancement factor (E s ) and the variable enhancement factor 
(£ v ); see text. 


a star and a habitable planet by up to 93%. By controlling 
for the natural extension due to larger eccentricity and 
obliquity, we found that their oscillations can extend the 
outer edge by up to 20% and never decrease it. Thus, the 
number of potentially habitable planets in the Galaxy may 
be larger than previously thought. 


We interpret our results to mean that planets with large and 
rapid obliquity oscillations are more likely to be habitable than 
those with negligible oscillations, such as Earth. This per- 
spective is at odds with the notion that the stability of Earth’s 
obliquity is important to the development of life. While it still 
may be true that rapid oscillations can be detrimental, and 



System Number 


FIG. 11 . A visualization for the HZ enhancement factors from Table 4. The height of the bars is the HZ enhancement factor, 
E s , for the complete simulations. The green box shows the fraction of that enhancement due to the variability of the system, E v . 
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FIG. 12. A comparison of how the eccentricity and obliquity impact the calculated outer edge of the HZ as determined by 
the temperature habitability index for the baseline opacity of 0.095. The top panels show the eccentricity (left) and obliquity 
(right) for the variable runs. The bottom panels show the same for the static cases. The error bars for the simulations 
represent the standard deviation in the eccentricity and the obliquity. Note there is little correlation with eccentricity but a 
strong correlation between obliquity and the outer edge of the HZ. When the variability in the runs is removed, the outer 
edge moves inward in all but four cases (see Table 4). Values for the linear fits are given in Table 5. (Color graphics 
available online at www.liebertonline.com/ast) 


certainly at some point obliquity cycles could be too large and 
rapid, our results clearly show that rapid obliquity evolution 
can be a boon for habitability. At the least, one should not rule 
out life on planets with rapid obliquity cycles. 

Our results are important for future telescopic searches 
for life, such as the Terrestrial Planet Finder (TPF). Al- 
though a final design has yet to be selected, TPF’ s mission is 
to directly image potentially habitable planets. These ob- 
servations depend critically on large star-planet separations 
in order to disentangle stellar light from reflected planetary 
light. Our results show that potentially habitable planets 
can exist at larger star-planet separations than previously 
appreciated, improving the odds that TPF can discover an 
inhabited planet. 


Table 5. Parameter Fits for the Linear Models 
in Figure 12 



Slope 

Intercept 

x 2 

Variable eccentricity 

0.16 

-0.19 

1.7 xlO 5 

Average eccentricity 

0.17 

-0.18 

3.5xl0 8 

Variable obliquity 

46 

-43 

4.1 

Average obliquity 

53 

-48 

91 


The calculation of the goodness of fit parameter, y 1 , along with 
the associated uncertainties, is described in the text. 


By necessity, our study was based on hypothetical planets. 
While our model systems are extreme from a Solar System 
point of view, they are relatively tame by exoplanet stan- 
dards. Eccentricities larger than 0.9 have been discovered 
(Naef et al, 2001; Jones et al, 2006; Tamuz et al, 2008), 
very large mutual inclinations are implied by the misalign- 
ment between some stars’ rotation axis and the orbital planet 
of a companion planet (Triaud et al, 2010; Naoz et al, 201 1), 
and eccentricity-inclination coupling can drive large and 
rapid oscillations, if both properties are large (Kozai, 1962; 
Takeda et al, 2007; Barnes et al, 2011). Thus, our results 
should not be viewed as an extreme possibility but rather as in 
the middle of a spectrum of spin-orbit coupling. 

Our approach has been simplified in several important 
ways. While our A-body simulations are the best represen- 
tation of possible orbits, our rotational model is simplified. 
A better approach would be to calculate the direct torques 
from all the bodies in the system and adjust the angular 
momentum distribution accordingly. Such a model is much 
more computationally expensive but could be incorporated 
into an A-body model without too much extra computational 
cost. Our EBM is also highly idealized, and future im- 
provements could include ocean/land dichotomies, the 
physics of glacier advancement and retreat, cloud physics, 
and ultimately even a three-dimensional global circulation 
model. Each of these additions, however, adds free 
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parameters to a model with very few constraints. While 
these features will improve realism, we will continue to 
suffer from a dearth of observational constraints. None- 
theless, as we move toward identifying planets worthy of 
detailed spectroscopic follow-up, such modeling could 
provide additional insight for prioritization. 

While our results demonstrate that planetary system ar- 
chitecture can influence the position of the HZ, it remains 
unclear how robust this influence is. For example, our 
planets all began with a spin rate of 24 h and an obliquity of 
23.5 degrees. How do different choices change the picture 
presented here? Future work should explore a range of ini- 
tial conditions and determine whether certain architecture 
always drives the planet into a particular obliquity cycle. If 
true, then we may be able to characterize a planet’s obliq- 
uity without direct measurements. While such a study was 
beyond the scope of this study, the possibility of tightly 
constraining obliquity is tantalizing and certainly worthy of 
a follow-up investigation. 

Our study suggests that rapid changes in obliquity and 
eccentricity increase the outer edge of the HZ. We quantify 
that relationship with linear trends in the enhancement 
factor with obliquity, but we did not find a threshold to 
achieve a specific quality that permits significant expansion. 
We blame the small number of systems we studied for this 
ambiguity and leave its identification for future work. We 
note that prior to running a simulation, it is very difficult to 
know how the orbital and rotational angular momenta will 
evolve; thus it could take considerable effort to produce a 
suite of architectures that suitably cover parameter space. 

Our study has shown how orbital architecture is a crucial 
factor when assessing planetary habitability. While previous 
work has mostly focused on static planetary properties, 
planets are expected to lie in multiplanet systems; hence the 
sequence of states must be considered. For the foreseeable 
future, we will have very few constraints on the properties 
of potentially habitable planets, and we therefore must 
leverage any information we have. 
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